function [Ss, rs] = find_spillovers_mech(Ss, Hs, rs, Prob, tol,max_It)
    global bet0 bet1 tau a b w e pi_theta N_a N_theta N_w N_e w_max w_min alpha myslope...
        xi xi1 ub lb pi_e H w theta Int_pi_w sigmae elasts Pr log_e_mu log_e_sigma phi sigma gamma rts r_bar fp r_base w_mult w_bar...
        
   
    pi_w = sum(Prob, 2);
    w_bar = dot(w, pi_w ); 
    
    % Spectral Algorithm to find market clearing prices
    x_0 = rs(1:end);
    x_1 = x_0 + (size(rs,1):-1:1)'/100;  

    f_0 = find_rents_dyn(x_0, Ss,  Hs, Prob);
    f_1 = find_rents_dyn(x_1, Ss,  Hs, Prob);

    max_iter = 100;
    for i = 1:max_iter 
        s = x_1(1:end) - x_0(1:end);
        y = f_1(1:end) - f_0(1:end);

        x_1 = x_0;
        f_1 = f_0;

        lambda = dot(s,y)/dot(y,y);

        step_size = 1;
        for iter_sub = 1:20

            x_0_t = x_1 - lambda * f_1 * step_size;

            if sum(x_0_t(1:end-1) > 0 ) == (length(x_0_t)-1)
                break
            else
                step_size = step_size * 0.9;
            end
        end
        x_0 = x_0_t;
        f_0 = find_rents_dyn(x_0, Ss,  Hs, Prob);


        diff = sqrt(dot(x_0 - x_1, x_0 - x_1));
        diff_f = sum(abs(f_0));

        if (diff < 1e-6) & (diff_f < 1e-3)
            rs = x_0;
            break
        elseif i == max_iter

            rs = [w_bar*0.7, w_bar*0.5, w_bar*0.2]';
            x_0 = rs(1:end);
            minf = @(x) sum((find_rents_dyn(x, Ss,  Hs, Prob)).^2);
            options = optimset('TolFun',1e-6,'TolX',1e-3);
            [minx, fval, exitflag] = fminsearch(minf, x_0(1:end), options );

            rs = minx; 
            if (fval > 1e-3 | isnan(fval) | (exitflag ~= 1)) & (it > 5 )
                fval, minx
                error("Did not find root")
            end
        end                
    end
    
end